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Realization of conically linear dispersion 
(termed Dirac cone as in Fig. [l](a)) has re- 
cently opened up exciting opportunities for high- 
performance devices that make use of the peculiar 
transport properties [Ij-^Gj of the massless carri- 
ers. A good example of current fashion is the 
heavily studied graphene, a single atomic layered 
graphite. It not only offers a prototype of Dirac 
physics in the field of condensed matter and ma- 
terials science [7J, but also provides a playground 
of various exotic phenomena [9-- ;12j . In the mean- 
time, numerous routes have been attempted to 
search for the next "graphene" [14— 20 J. Despite 
these efforts, to date there is still no simple guide- 
line to predict and engineer such massless parti- 
cles in materials. Here, we propose a theoret- 
ical recipe to create Dirac cones into anyone's 
favorite materials. The method allows to tailor 
the properties, such as anisotropy and quantity, 
in any effective one-band two-dimensional lattice. 
We demonstrate the validity of our theory with 
two examples on the square lattice, an "unlikely" 
candidate hosting Dirac cones, and show that a 
graphene-like low-energy electronic structure can 
be realized. The proposed recipe can be applied 
in real materials via introduction of vacancy, sub- 
stitution or intercalation, and also extended to 
photonic crystal [21j, molecular array [22j, and 
cold atoms systems [23j. 

Let's first review an interesting observation of the 
Dirac cones [24 in graphene via angle-resolved photoe- 
mission spectroscopy (ARPES) [25ti27 . since the resolu- 
tion of it would lead naturally to our proposed method. 
Figure [TJb) shows an representative ARPES observation 
of the Dirac cone using out-of-plane polarized light. In- 
triguingly, the observed cone appears incomplete even 
though a standard theory would indicate a complete 
Dirac cone, for example given by the intense red bands 
in Fig.jljd). This vanishing intensity is typically consid- 
ered a "matrix element effect" of the measurement, and 
indicates a perfect destructive quantum mechanical in- 
terference [28 between the two carbon atoms in the unit 
cell of the honeycomb lattice shown in Fig. [ijc). 

An alternative and more straightforward perspective 
is to consider the honeycomb lattice as a triangular lat- 



tice with a periodic vacancy, as shown in Fig.jTJe). Since 
there is only one (or sometimes zero) atom in this unit 
cell, the matrix element would just be the simple atomic 
form factor and the vanishing intensity is thus explicitly 
incorporated in the corresponding one-particle spectral 
function. Indeed, Fig. [ijf-h) shows the resulting "un- 
folded" spectral function [29j with a "incomplete" Dirac 
cone having vanishing spectral weight in the lower part 
of the cone along the KgMn path, and in the upper part 
along the FKg path, in perfect agreement with the exper- 
imental observation. 

This alternative "one-carbon" picture offers an inter- 
esting new way to understand the electronic structure of 
graphene, particularly the formation of the Dirac cones. 
Fig. [2] illustrates this with a reduced Hamiltonian, for 
clarity, that covers only the C-pz orbitals that define the 
relevant red band in Fig. [l] Starting with a triangu- 
lar lattice of carbon atoms with the same nearest inter- 
atomic coupling t as in graphene, the corresponding band 
structure consists of a simple dispersion [Fig. |2|a)] and 
an almost circular Fermi surface [Fig. [2|d)] when the or- 
bital is half-filled. Upon raising the energy of the ordered 
vacancy sites by £, the system is driven into a charge 
density wave (CDW) state, with most part of the Fermi 
surfaces (b)&(c) gapped out, leaving six Dirac cones in 
the original one-carbon Brillouin zone (BZ) [denoted by 
solid black boundary lines in (d)(e)&(f)]. As s grows, the 
effects is further enhanced and the cone becomes more 
and more symmetric [Fig. [2|j)(k)]. Finally, the perfectly 
symmetric cone of graphene is reproduced as the vacancy 
sites become forbidden, e ^ 00. 

Deeper physical and mathematical insights can be 
obtained by investigating the analytic structure of the 
above CDW potential in the momentum space: 
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where r and denotes the atomic sites, k the crystal 
momentum and q the CDW wave vector. With only a 
local energy increase at the vacancy sites Vry = eSry , 
Vk^k-\-q displays a very specific form within the subspace 
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of the three states that couple together: 




yk,k-\-q 



FIG. 1: New perspective of graphene lattice. (a) 

Schematic diagram of linear dispersion of Dirac cones in 
graphene at the Brillouin zone corner (K point). On the base, 
the red circle represents the energy contour at 0.4 eV above 
the cone Fermi energy, (b) ARPES data from Ref. [27^ at 0.4 
eV above Dirac point, (c) The graphene honeycomb lattice 
of the blue unit cell (as the super cell), (d) The first-principle 
band structure with color code: red for pz and green for px 
& Py orbit als. The zero energy is chosen at the Fermi energy 
of graphene. The bands of two different color types are fully 
decoupled because the coupling is forbidden by the mirror 
symmetry with respect to the graphene sheet, (e) A refer- 
ence triangular lattice of the black unit cell (as the normal 
cell) with the imaginary gray carbon, located at the hexagon 
center and treated as an vacancy, (f) The first-principle un- 
folded band structure in the normal cell basis. The thickness 
represent the magnitude of spectral weight of Green's func- 
tion. The subscripts "s" and "n" are used to distinguish the 
high symmetry points defined for the super and normal cells, 
(g) and (h) represent the unfolded energy iso-surfaces at en- 
ergy = ±1 eV. 




diagonalization 




Applying such 
Fig. [2];a)(d)&(g) 



(2) 

I coupling to three states in 
that happen to be degenerate, 
one finds a two-fold degeneracy is preserved in the 
resulting eigenvalues, forming the Dirac points, the tips 
of the Dirac cones. Furthermore, with a slight deviation 
dk from these k points, the energy of the coupled states 
in (a) would differs by small amounts a and 6, both 
proportional to 6k. 
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Taking \a\ < \b\ <C e for convenience, one finds linear 
dispersion developing around the Dirac point: f + 
— Iff}. Within such a one-carbon picture, we 
have demonstrated the formation of Dirac cones through 
an induced CDW of specifically structured CDW poten- 
tial. The above analysis also locates the Dirac cones in 
momentum space from the reference CDW-free system 
before the CDW q vectors to be applied are given. 

The vanishing spectral intensity in Fig. [l] is also easily 
understood with a similar analysis. Along the k path on 
which two of the three coupled states remain degenerate 
(a = 0), for example two /c-points on the KnKgMg path 
[gray dotted lines in Fig. |2|b)(e)] coupled to the third on 
the FKgMn path [gray dashed line in Fig. |2|b)(e)], one 
of the resulting eigenvalues remains unchanged (zero): 
{3 + f,^,0}, and the corresponding eigenvector con- 
sists of anti-bonding superposition of only the two orig- 
inally degenerate states. Consequently, this band will 
not be folded to the third k point (on the gray dashed 
line), where the corresponding spectral intensity must 
then vanish in Fig{TJg)&(h). The absence of folding in- 
tensity may be more clearly visualized in Fig. [2|a)-(c), 
where the band along the red path in (a) never appears 
in the green path, against naive expectation. In essence, 
within this one-carbon picture, the vanishing spectral in- 
tensity is intimately tied to the formation of the Dirac 
cones. 

The above generic description suggests a powerful and 
general recipe to create Dirac cones in any quasi-2D ma- 
terials. It would allow adding the exotic physics of Dirac 
particles and enhanced transport property into one's fa- 
vorite materials that already hosts other useful character- 
istics. This is achieved by inducing a CDW state through 
introduction of impurities like vacancy, substitution, or 
intercalation. These kinds of impurities generates im- 
purity potential dominant ly local at the impurity sites, 
thus giving a nearly constant CDW coupling in momen- 
tum space similar to Eq. |2] The sample should then be 



3 



a 8/t = b 8/t = 4 c 8/t=10 




FIG. 2: Modeling of Dirac cone creation induced by CDW effects on the triangular lattice, (a)-(c), unfolded band 
structure for e/t — 0, 4, 10^. (d) - (f) represent energy isosurfaces in the normal cell basis, namely Afe,fc(<x;), for three different 
e/t. Different color stands for different energy cut with range from -0.8 t to 0.4 t. In order to show details of non-uniform 
spectral weight distribution, the yellow dashed square region in (e) is enlarged into an inset between the second and third row. 
(g) - (i) are the energy iso-surfaces plotted in super cell Brillouin zone. To directly visualize Dirac cones, (j) and (k) are the 
energy dispersion of the cyan and black dashed square in (h) and (i) respectively. The very flat band on top cone in (j) is 
responsible for the hexagonal shape of spectral function in (e) and (h). In (b) and (e), the two degenerate gray dotted k path 
is coupled to the third gray dashed path by CDW potential. The red dashed curve in (b) refers to the missing folded spectral 
due to the special format of the potential. 



synthesized/annealed to homogeneous distribution of the 
impurity, such that at higher enough concentration the 
impurities are expected to order reasonably well. The key 
CDW wave vector would then be controlled by the impu- 
rity concentration, to couple M (M > 3) /c-points in the 
clean system. Since geometrically it is always possible 
to couple at least three degenerate states with a regular 
dispersion, the above derivation would dictate the ap- 
pearance of Dirac cones around the M /c-points of these 
degenerate states. In principles there can be multiple 
sets of such coupled /c-points, in which case the number 
of Dirac cones would multiply. As long as the structure 
of the materials is stable against the not-so-small amount 
of impurity, this recipe offers a controllable and robust 
means to generate massless Dirac particles in almost any 
materials of interest. 

To demonstrate the general validity of this recipe, let 



us take a 2D one-band system with a common square lat- 
tice and introduce 33.3% substitution of the atoms uni- 
formly across the system. Figure [3] demonstrates creation 
of six anisotropic Dirac cones in the original BZ (two in 
the new reduced BZ). In the case with weaker impurity 
potential e (Fig. [3jc)&(f)), three regular "electron pock- 
ets" remain at the chemical potential. In general, contri- 
bution to transport properties from the normal massive 
carriers on these pockets should be overwhelmed by those 
of the massless Dirac carriers, and thus does not cause 
any serious concern. These normal carriers can often be 
removed by gapping out the pockets with a stronger im- 
purity potential (for example with vacancy), as shown 
in Fig. [3]^d)&(g). Obviously, the closer the reduced BZ 
is to the original Fermi surface, the more effectively the 
impurity can gap out these Fermi pockets. 

One very interesting utilization of our recipe is to try to 
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8/t = 8/t = 5 s/t = 20 




00 - 0.58 t 
03 = 0.5 t 
03 = 0.4 t 
03 = 0. 32 t 
03 = 0.24 t 



FIG. 3: Dirac cones on square lattice of 33% ordered impurities (M=3). (a) The original square lattice of black 
normal cell and the blue super cell of {A = (1,1), B = (1,2)}. The zero energy and next nearest hopping are chosen as the 
on-site energy of the unperturbed orange atoms and 0.4 t respectively. The on-site energy of the gray atom is lifted to e. (b) - 
(d) are the unfolded energy isosurfaces plotted in normal cell basis for s/t =0, 5, 20. Different color represents different energy 
cut with range from 0.24 t to 0.58 t. (e) - (g) are energy isosurface in super cell basis, (h) is the energy dispersion of the region 
of cyan dashed square in (g). 



s/t = 



8/t=5 



s/t = 10 




03 = -0.82 t 
03 = -0.92 t 
co = -1.0 t 

03 = -1.12t 
03 = -1.2 t 



FIG. 4: Dirac cones on square lattice of 25% ordered impurities (M=4). (a) The supercell of {A = (2, 1), 5 = (0, 2)}. 

The graph convention follows Fig. reffig:fig3. (b) - (d) are the energy isosurfaces plotted in normal cell basis for s/t =0, 5, 20. 
Different color represents different energy cut with range from -1.2 t to -0.82 t. (e) - (g) are energy isosurfaces in super cell 
basis, (h) is the energy dispersion of the region of cyan dashed square in (g). 



create a graphene-like electronic structure from a square 
lattice, rather than hexagonal structure. Figure [4] shows 
that introducing 25% uniformly distributed impurities 
can indeed create six Dirac cones in the original BZ that 
resemble very much those in graphene. Out of four cou- 
pled k points, only three of them are degenerate (M=3). 
The state at the fourth /c-point has a much higher energy 
and thus does not affect the other three in any significant 



manner. 

Finally, we clarify a few common practical aspects in 
real materials. First, it is obviously impossible to have a 
perfect order in the impurity. Luckily, the linearly van- 
ishing density of states to the Dirac point and its massless 
feature would render the disorder effects of weak imper- 
fection rather insignificant. In essence, Dirac carriers are 
quite immune to disorder effects in general [8 . Second, 
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our analysis above takes a simple limit that the impurity 
potential is dominantly local in the impurity site. While 
this is indeed the case for almost all the impurities (and 
certainly can be controlled by choosing the atoms for sub- 
stitution), a small non-local influence of the impurity is 
expected in real materials. Such a small non-local effect 
would introduce a small /c-dependence of the coupling, 
and thus a small variation of the off-diagonal elements 
of Eq. |2] This in turn would introduce a second order 
correction S to the resulting eigen-energy and possibly 
lift the degeneracy of the Dirac point with a small gap. 
At the temperature/energy scale larger than S one would 
still observe Dirac linear dispersion. On the other hand, 
this offers an additional engineering degree of freedom to 
create a controlled semi-conducting gap that itself can 
have great deal of applications [8 . 

Third, for system with more 3D like dispersion, our 
mathematical argument would produce Dirac cone struc- 
ture in the 2D CDW plane, but with Dirac point con- 
nected in the direction perpendicular to the plane. This 
leads to an interesting situation of carriers having zero 
mass in the plane, but an almost infinite mass in the out- 
of-plane direction, essentially promoting a 2D character- 
istic. Finally, the created Dirac cone might not always 
occurs at the chemical potential in real material. This 
can be overcome by the standard technique of doping or 
applying gate voltage. 
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Supplementary information 



In the supplementary information, we will provide computational details about the ffist-principle results graphene 
and band structure unfolding mentioned in the main text. Standard density functional theory (DFT) calculation 
is performed by using full potential linearized augmented plane wave method with local density approximation im- 
plemented in the WIEN2k package [1 version 10.1. To simulate the quasi- two-dimensional properties of graphene 
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of honeycomb lattice, the unit ceh, containing two atoms, is chosen with 1.42 A inter-carbon distance and 10 A 
inter-layer distance vertically. This forms a crystal structure of P6/mmm space group. All input settings follow the 
default values as RmtKmax = 7 and Imax = 10. The k point mesh of 21 x 21 x 4 is used to reach convergent ground 
state density. Then the symmetry-respecting Wannier functions of carbon s and p characters would be constructed 
within [-20:30] eV low energy Hilbert space based on first-principles dispersion [2 . 

As mentioned in the main text, it is physically meaningful to represent the electronic band structure in the reference 
system corresponding to the one-carbon unit cell on triangular lattice. Indeed, this unfolding technique helps us to 
reveal the important symmetry breaking effect: the vacancy-induced charge density wave leads to the formation of 
Dirac cones. The basic idea of our unfolding method [3] is to simply represent the energy {uj) dependent one-particle 
spectral function, imaginary part of Green's function, of the real system (honeycomb graphene) by using the basis 
from reference system (triangular carbon lattice): Akn,kn{^) = \ {^^\^J)\^^kj,kj{^)^ where K /k denotes the 

crystal momentum of the graphene/triangular lattice, J the band index, and n the Wannier orbital index (like s and 
p orbitals). This change of basis is made simple with the use of the above symmetry-respecting Wannier functions 
[2 . For the demonstration purpose in the main text, we show both Akj^kj{^) and Akn,kn{^) to illustrate the clearer 
morphology and spectral weight distribution of the Dirac cone respectively. 
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